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Introduction 

While there has been much recent progress in simulating nonlinear aeroelastic systems, and 
in predicting many of the aeroelastic phenomena of concern in transport aircraft design (i.e. 
transonic flutter buckets), the utility of a simulation in generating an understanding of the 
flutter behavior is limited. This is due in part to the high cost of generating these 
simulations, and the implied limitation on the number of conditions that can be analyzed, 
but there are also some difficulties introduced by the very nature of a simulation. Flutter 
engineers have traditionally worked in the frequency domain, and are accustomed to 
describing the flutter behavior of an airplane in terms of its V-G and V-F (or Q-G and Q-F) 
plots and flutter mode shapes. While the V-G and V-F plots give information about how 
the dynamic response of an airplane changes as the airspeed is increased, the simulation 
only gives information about one isolated condition (Mach, airspeed, altitude, etc.). 
Therefore, where a traditional flutter analysis can let the engineer determine an airspeed at 
which an airplane becomes unstable, while a simulation only serves as a “binary check:” 
either the airplane is fluttering at this condition, or it is not. 

In this document, a new technique is described in which system identification is used to 
easily extract modal frequencies and damping ratios from simulation time histories, and 
shows how the identified parameters can be used to determine the variation in frequency 
and damping ratio as the airspeed is changed. This technique not only provides the flutter 
engineer with added insight into the aeroelastic behavior of the airplane, but it allows 
calculation of flutter mode shapes, and allows estimation of flutter boundaries while 
minimizing the number of simulations required. 


Aeroelastic Equations of Motion 

For the purposes of this document, it is useful to consider the aeroelastic system as a 
separate structural dynamic system (which depends only on the structure’s stiffness and 
mass properties) and aerodynamic system (which depends on the flow parameters). These 
systems are then conceptually coupled into an aeroelastic system using a feedback loop, 
where the generalized deflections from the structural system drive the aerodynamic model, 
and the generalized forces from the aerodynamic model in turn drive the structural model. 
This is shown schematically in Figure 1. 



Generalized 

Deflections 


Figure 1: Block Diagram Representation of the Aeroelastic System. Note 

that Inputs and Outputs Are Vector Quantities. 


Before considering the coupled aeroelastic system, it is useful to define the form of the 
dynamic equations governing the structural dynamic system and the aerodynamic system 
separately. Invoking a small perturbation assumption, the dynamics of each system can be 
approximated using a linear first-order matrix equation. The resulting state-space equations 
for the structural model are given below. 


*,VI = A s*i + B sfi 
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where Xj is a vector containing the (structural) generalized deflection and generalized 
velocity information describing the airplane motion, and is the generalized aerodynamic 
force. 

The first order differential equation governing the aerodynamic system can be written as 


%i+l A a <^/ + B z Xj 

f i =qC£ i +qDx i 


where ^ is a vector containing a description of the state of the aerodynamic system. Note 
that this state vector is in general extremely large, containing a degree of freedom for every 
degree of freedom of the CFD model. For example, in a potential based code, would 
contain the potentials for every grid point in the CFD model. 

Combining These Equations Leads to 

* 1+1 = ( A S + qB s D)x t + qB s C4i 

6+1 = A aZi +B a X i 


or 



~A S + qB s D qB s C!M 

. B a A a Mi 


( 4 ) 


This equation representing the coupled aeroelastic system will form the basis of the flutter 
boundary identification algorithm described below. 


Algorithm 

There are two significant features of the method presented here. The first is the application 
of system identification methods to efficiently and accurately estimate the dynamic behavior 
(including modal frequencies and damping ratios) of the coupled aeroelastic system. The 
second contribution is a technique for expanding the resulting information (which is 
applicable only to the exact combination of Mach number and dynamic pressure simulated 
in the CFD code) into a description of how the system dynamics vary with changes in 
dynamic pressure. 

Identification of System Dynamics 

Consider the dynamic system described in equation (4) above. Note that this representation 
appends the aerodynamic states to the structural states to form the coupled aeroelastic 
system. Strictly speaking, the number of aerodynamic states is extremely high, and scales 
with the number of grid points in the CFD model. In practice, however, the dynamics of 
the aerodynamic system can be adequately represented by a much smaller set of degrees of 
freedom that contain essentially the same information (conceptually similar to aerodynamic 
lags). The number of “important” states is often very small, and some very important 
information about the aeroelastic system can be obtained by treating the aerodynamic states 
as second order effects and simply neglecting them. Truncating the aerodynamic states 
results in an approximate dynamic equation of the form 

*/+i ={A s +qB s D)x i (5) 
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Note that a large part of the aerodynamic force is accounted for in a quasi-unsteady fashion 
in the term qB s D, while the increment to fully unsteady aerodynamics is neglected. This 
can be written as 


* i+i = ^ae(^) x i ( 6 ) 

where A^q) is a state transition matrix for the aeroelastic system, and is obviously 
dependent on dynamic pressure. The problem then is to estimate the matrix A^q) given 
the time history of the generalized deflections x ; . 

This is accomplished using a simple least squares approach. First, the history of the 
generalized displacements and generalized velocities are assembled into a matrix: 

-*?=[*« x »-i] 

K = hi - *„] (7) 

and the dynamic equation of the aeroelastic system over the time history can be written in 
matrix form as 

K = A AE (.q)X° n ( 8 ) 

Since the aerodynamic lag terms are neglected, and since the system is not perfectly linear, 
there is not an exact solution for A AE (q), but it can be estimated in a least squares sense: 

A AE [q) = X'X T {x°X T ) 

This can also be written as 

A A £(?) = £/„W 

where the matrices U n and V n are defined as 

U n = Z*; + \ X J 

i = o 

K = 

i = 0 

In this approach, the state transition matrix A^ E (q) can be estimated at every iteration as the 
simulation progresses, and by tracking the eigenvalues (frequency and damping ratio) of 
the state transition matrix, an excellent feeling for the convergence of the solution can be 
obtained. This makes it easy to determine when to stop CFD solutions, so CPU time is not 
wasted. 

Interpolation/Extrapolation to Different Conditions 

If a the state transition matrix at is available at at least two dynamic pressures, a flutter 
analysis can be performed by using linear interpolation/ extrapolation on a term-by-term 
basis to estimate a state transition matrix at any dynamic pressure of interest: 


( 9 ) 

( 10 ) 

(H) 
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( 12 ) 


A ae (?) S - 22 — 1 A ae (?, ) + 3-21- Aae ( ?2 ) 

#2 #1 *?2 #1 

This matrix can be inexpensively estimated at a large number of dynamic pressures, and its 
eigenvalues can be calculated. These can be transformed into the continuous-time domain 
and be used to calculate modal frequencies and damping ratios of the aeroelastic system. 
This information can be tabulated for all the dynamic pressures, and a set of Q-G and Q-F 
plots can be constructed, allowing flutter interactions to be identified, and flutter crossings 
to be estimated. It should be emphasized that the resulting plots are exact only at q=q, and 
q=q,, and are “unmatched” approximations at other dynamic pressures. 

In this report, a single aeroelastic CFD analysis was performed at q 2 , and q, was taken to 
be zero. The known frequencies and damping ratios of the structural dynamic system were 
then used to construct a state transition matrix for the wind-off condition A ae ( 0). In 
general, to avoid excessive interpolation/extrapolation errors, multiple aeroelastic analyses 
might be required, and piecewise linear or higher order polynomial interpolation might be 
used. A preliminaiy analysis using the wind-off dynamics and a state transition matrix at 
an initial guess of dynamic pressure could also be used to guide the engineer’s choice of 
additional dynamic pressures to run to ensure adequate accuracy to predict the flutter 
crossing. 


Results: FSM Analysis 

The algorithms described above have been applied to several configurations, including the 
AGARD 445.6 wing and the HSR Flexible Semispan Model (FSM). A typical set of 
results from the FSM configuration are shown here for illustrative purposes. 

The structural model used in this example is the Long Beach version of the correlated Ml 8- 
5 Finite element model. This is documented thoroughly in Reference 1, and a plan view of 
the FEM is shown in Figure 1. In order to perform the aeroelastic solutions, the vibration 
modes of the Finite Element Model are calculated, and used as input into the Boeing Long 
Beach version of CFL3D-AE (CFL3D-AE-BA), which is described in Reference 2. The 
Euler grid used in this example is shown in Figure 2. 

In order to compute a flutter simulation, a steady-state aeroelastic solution must first be 
obtained in order to accurately capture the mean steady flow properties. This is 

accomplished in two steps. First, a rigid steady-state solution is computed, and then the 
solution is restarted with the effects of structural flexibility included. The convergence 
plots of the static aeroelastic solution (in terms of modal generalized displacements of the 
dominant modes) are shown in Figure 3. The flutter analysis is then performed by 
restarting the aeroelastic solution, and by imparting a nonzero generalized velocity to each 
of the modes. In this case, the analysis was performed at Mach 0.8, Q=250 PSF, and 
angle of attack of 2 degrees. Ten modes were used in the aeroelastic solution. The 
resulting time histories of generalized deflections then represent the impulse response of the 
coupled aeroelastic system, and are shown in Figure 4. Only modes 1, 3, 4, and 6 are 
shown because these show the most interesting dynamics. 

On inspecting Figure 4, it can clearly be seen that we are very close to neutrally stable in at 
least one aeroelastic mode (the generalized deflections of Modes 1 and 3 are virtually 
undamped). There is also some interaction in Modes 4 and 6, since the presence of 
different frequency components is evident in the time histories. This is expected since the 
generalized deflections represent the natural modes of the wind-off system, and once the 
dynamic pressure is raised above zero, coupling between generalized coordinates 
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introduces multiple frequency components. All that can really be ascertained from these 
plots, however, is that this condition is not significantly unstable. 

Figures 5 and 6 show the application of the system identification technique for estimating 
the eigenvalues of the aeroelastic system. Figure 5 shows the convergence of the modal 
frequency estimates, and Figure 6 shows the convergence of the modal damping estimates. 
The algorithm exhibits some chaotic behavior for about the first 1 00 iterations, but has no 
trouble converging on frequency and damping estimates in about 200-250 iterations. This 
represents about 5 cycles of the dominant aeroelastic mode. The results of this analysis are 
consistent with our visual interpretation of Figure 4: there is a marginally stable mode at a 
frequency of about 13-14 Hz, and a slightly more stable mode at a frequency of about 18 
Hz. 

The FSM results using CFD (CFL3D-AE-BA) and system identification are expanded into 
Q-G and Q-F plots in Figures 7 and 9, and the results of a linear flutter analysis are 
presented for comparison in Figures 8 and 10. From Figure 7, it can clearly be seen that 
the first aeroelastic mode rises in frequency as dynamic pressure is increased, and crosses 
the third aeroelastic mode at a dynamic pressure of about 200 PSF and the fourth 
aeroelastic mode at a dynamic pressure of 300 PSF. Above 300 PSF, this mode’s 
frequency starts dropping off. The other significant effect is that the sixth aeroelastic 
mode’s frequency drops as the Q is increased, and the frequency shows a distinct inflection 
at about 300 PSF. 

Compare the frequency behavior from the nonlinear solution with the results of the linear 
flutter analysis in Figure 8. Again, we see that the frequency of the first mode increases 
with Q, and that the first and third modes cross at roughly 180 PSF. The first and fourth 
modes interact at about 275-325 PSF, after which the frequency associated with the first 
mode begins to fall off. We also see that the sixth elastic modal frequency drops as Q is 
increased, and that there is a distinct inflection at about 320 PSF. The comparison between 
the Q-F curves deteriorates above about 400 PSF. 

Now consider the damping behavior as dynamic pressure is increased. In the nonlinear 
results (Figure 9), the existence of three aeroelastic mechanisms can be seen. The first is a 
hump mode that peaks at a damping value just below zero at about 250 PSF. By referring 
to the Q-F plot, it can be seen that this mode corresponds to the third elastic mode, and has 
a frequency of about 13 Hz. Figure 10 shows a similar mechanism for the linear analysis, 
but the peak damping value is about 5% beyond the stability boundary. 

The second aeroelastic mechanism seen in the nonlinear results is a hard flutter crossing 
occurring at around 260 PSF. This mode corresponds to the fourth elastic mode, and has a 
frequency of about 1 8 Hz. Inspecting the linear analysis shows a similar flutter mechanism 
with a crossing at about 280 PSF. 

The third mechanism in the nonlinear results is again a hump mode that peaks at nearly 500 
PSF, significantly below the stability boundary. Since the results at Q=500 PSF are 
extrapolated from Q=250 PSF, there is a low degree of confidence in this mechanism. The 
linear results show action in this mechanism, but the interaction between mechanisms is 
significantly different at the high dynamic pressures involved. 


Conclusions 

The techniques presented in this document have been used to estimate the flutter boundaries 
of the FSM configuration, and have provided some valuable insight into the aeroelastic 
interatcions involved. This approach offers three main advantages over other techniques: 

1. It gives the engineer insight into the adequacy of the simulation for identifying system 
dynamics (through convergence plots of modal frequencies and dampings). 
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2. It allows interpolation/extrapolation of the identified dynamics to predict the flutter 
behavior at other dynamic pressures, and to see the aeroelastic interactions leading to 
flutter modes. 

3. It allows the potential for calculating flutter mode shapes from the eigenvectors of the 
identified state transition matrix. 

Some caveats on these results are also in order. To date, the validation of these techniques 
is fairly sparse, and questions about the range of applicability of interpolated/extrapolated 
solutions and the application to highly nonlinear conditions (i.e. Mach 0.95) still need to be 
answered. Further validation will take place as the nonlinear analysis of the FSM 
configuration continues. 

The source code of the program that performs this analysis is shown in the Appendix. 
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Appendix: Source Code 

program flutter 
c maximum number of modes 
parameter (nmax=20) 
implicit real*8 (a-h,p-z) 

real *8 d(nmax) ,v(nmax) ,f (nmax) ,dt,u0,cbar, stm(nmax, nmax) 
real*8 freq (nmax) , damp (nmax) ,wr (nmax) ,wi (nmax) 
complex* 16 e 
character* 132 f n, line 

ff = 1.000 
nq = 300 
qmaxq = 3.0 

c CAP-TSD style MCFOUT file containing modal time history data 
iin = 9 

c output file for convergence of modal frequencies and damping ratios, 
iconv =10 

open ( iconv , f i 1 e= ' f reqdamp . conv 1 ) 
c flutter crossing file 
iflut = 11 

write(*,*) 'enter the filename of the CAP-TSD style MCFOUT file.' 
read(*, 1 (al32) ' ,end=700) fn 
write(*,*) 'opening file:' 
write (*,' (al32) ' ) fn 

open (iin, f ile=fn, status= 1 old 1 , err=999) 
read(iin2,*) nmmax 

read(iin2,*) itmax 

do i = 1,28 
read (iin, *) 
enddo 

read (iin, 1 (28x,fl2.6,53x,fll.6,llx,fl2.6) ') dt,q0,u0 
write (* , ' (28x,fl2.6,53x,fll.6,llx,fl2.6) ’) dt,q0,u0 
do i = 1,5 
read (iin, *) 
enddo 

read (iin, ' (59x, il3) ' ) nm 
write(* , ' (59x, il3 ) ' ) nm 
do i = 1,5 
read (iin, * ) 
enddo 

read(iin, ' (59x, f 13 .6) ' ) cbar 
write(* , ' (59x, f 13 .6) ' ) cbar 
dtau = dt*cbar/u0 
3 read (iin, ' (al32) ' ) line 

if ( index ( line, ’ gmass '). eq. 0) goto 3 
read (iin, *) 
read (iin, *) 
do i = 1 , nm 

read (iin, ' (27x, 2e20 . 12) ' ) damp ( i) , f req (i) 
write (*, ' (27x, 2e20 . 12) ' ) damp (i) , f req(i) 
enddo 

do i = l,2*nm 
do j = 1, 2*nm 
stm (i, j ) = 0 
enddo 
enddo 
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do i - l,nm 

aaa = - damp(i) *freq(i) 
bbb = sqrt (f req (i) **2 - aaa*aaa) 
ec = exp (aaa*dtau) *cos (bbb*dtau) 
es = exp (aaa*dtau) *sin (bbb*dtau) 
stm(2*i-l, 2*i-l) = ec - aaa*es/bbb 
stm(2*i-l,2*i ) = es/bbb 

stm(2*i ,2*i-l) = * (aaa*aaa + bbb *bbb )* es/bbb 
stm(2*i , 2*i ) = ec + aaa*es/bbb 

enddo 

write (*, *) ' stm: ' 
do i = l,2*nm 

write(* f ' (255el5 .5) ') (stm{i,j) , j=l,2*nm) 
enddo 

do i = 1,5 
read (iin, *) 
enddo 
iter = 0 
10 idone = 1 

read (iin, *,end=20) 
read (iin, *,end=20) 
read (iin, *,end=20) 
write(*,*) 'iter = ',iter 
do i = l,nm 

read (iin, ' (4x, 4 (3x, e21 . 12) ) ' , end=20) d(i) ,v(i) ,b, f (i) 
write (*, 1 (7x,e21.12,3x,e21.12) ') d(i) ,v(i) ,f(i) 
enddo 

write (25, ' (255el5.5) ') (d (i) , v (i) , f (i) , i=l, ran) 

iter = iter+1 
idone = 0 

if (iter.gt.itmax) idone = 1 
20 continue 

call sysid (iter, d, v, min (nmmax,nm) ,dtau, ff , stm,nmax, 

. iconv, qO, nq, qmaxq, if lut, idone) 

if (idone. eq.O) goto 10 

close ( iin) 
close (iconv) 
goto 1 

700 continue 
stop 

999 write (*,*) 'error.' 
stop 
end 

SUBROUTINE cholsl (a, n, np,m, p, b, x) 

INTEGER n,np 

REAL* 8 a(np,np) ,b(np,m) ,p(n) ,x(np,m) 

INTEGER i,k 
REAL *8 sum 
do 100 kk = l,m 
do 12 i=l,n 
sum=b (i,kk) 
do 11 k=i -1, 1,-1 

sum=sum-a (i , k) *x (k, kk) 

11 continue 
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x ( i , kk) =sum/p ( i ) 
continue 
do 14 i=n, 1,-1 
sum=x(i,kk> 
do 13 k=i+l,n 

sum=sum-a (k, i) *x (k, kk) 

13 continue 

x ( i , kk) =sum/p ( i ) 

14 continue 

100 continue 

return 

END 

C (C) Copr. 1986-92 Numerical Recipes Software 41.9;%{$1 
SUBROUTINE choldc (a, n, np,p, sing) 

INTEGER n, np 

REAL* 8 a (np,np) ,p (n) 

INTEGER i, j , k 
LOGICAL s ing 
REAL* 8 Sum 
sing = .false, 
do 13 i=l,n 
do 12 j=i,n 
sum=a (i, j ) 
do 11 k=i-l, 1,-1 

sum=sum-a(i,k) *a(j,k) 

11 continue 

if (i.eq. j) then 

if (sum. le.O . ) then 
sing = .true, 
return 
endif 

p (i) =sqrt (sum) 
else 

a { j , i) =sum/p(i) 
endif 

12 continue 

13 continue 
return 
END 

C (C) Copr. 1986-92 Numerical Recipes Software 41.9;%{$1 
SUBROUTINE elmhes (a, n, np) 
implicit real*8 (a-h,p-z) 

INTEGER n,np 
REAL* 8 a(np,np) 

INTEGER i , j , m 
REAL*8 x, y 
do 17 m=2,n-l 
x=0 . 
i=m 

do 11 j=m,n 

if (abs (a ( j ,m- 1) ) .gt.abs(x)) then 
x=a(j ,m-l) 
i=j 
endif 

11 continue 

if (i.ne.m) then 
do 12 j=m-l,n 
y=a ( i , j ) 
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a(i, j) =a(m, j) 
a (m, j ) =y 
continue 
do 13 j=l,n 
y=a(j,i) 
a ( j , i) =a ( j ,m) 
a(j ,m)=y 

13 continue 
endif 

if (x.ne.O.) then 
do 16 i=m+l,n 
y=a (i,m-l) 
if (y .ne. 0 . ) then 
y=y/x 
a (i,ra-l) =y 
do 14 j=m,n 

a(i,j)=a(i,j) -y*a(m,j) 

14 continue 

do 15 j=l,n 

a ( j ,m) =a ( j ,m) +y*a(j ,i) 

15 continue 
endif 

16 continue 
endif 

17 continue 
return 
END 

C (C) Copr. 1986-92 Numerical Recipes Software 41.9;%{$1. 
SUBROUTINE hqr ( a , n , np , wr , wi ) 
implicit real*8 (a-h,p-z) 

INTEGER n,np 

REAL*8 a (np, np) , wi (np) , wr (np) 

INTEGER i, its, j ,k, l,m, nn 
REAL*8 anorm,p,q, r, s, t,u, v,w,x,y, z 
anorm^abs (a (1, 1) ) 
do 12 i=2,n 
do 11 j=i-l,n 

anorm=anorm+abs ( a ( i , j ) ) 

11 continue 

12 continue 


nn=n 
t=0 . 

1 if (nn.ge. 1) then 

its=0 

2 do 13 l=nn, 2 , -1 

s=abs (a(l-l,l-l)) +abs (a (1 , 1) ) 
if (s .eq. 0 . ) s=anorm 
if (abs (a (1 , 1 - 1) ) +s .eq. s) goto 3 
13 continue 

1=1 

3 x=a(nn,nn) 

if (l.eq.nn) then 
wr (nn) =x+t 
wi (nn) =0 . 
nn=nn-l 
else 

y=a (nn- 1 , nn- 1) 
w=a(nn,nn-l) *a(nn-l,nn) 
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c 


14 


15 
4 

16 


if (1 .eq.nn- 1) then 
p=0 .5* (y-x) 
q=p**2+w 
z=sqrt (abs (q) ) 
x=x+t 

if (q.ge.O. ) then 
z=p+sign(z,p) 
wr (nn) =x+z 
wr (nn- 1) =wr (nn) 
if ( z . ne . 0 - ) wr (nn) =x- w/z 
wi (nn) =0 . 
wi (nn- 1) =0 . 
else 

wr (nn) =x+p 
wr (nn- 1) =wr (nn) 
wi (nn) =z 
wi (nn- 1) =- z 
endif 
nn=nn-2 
else 

if (its .eq. 30) pause 'too many iterations in hqr' 
if (its .eq. 1000) return 
if ( (its/10) *10 . eq. its) then 
t— t+x 

do 14 i=l , nn 

a (i, i) =a (i, i) -x 
continue 

s=abs (a(nn,nn-l) )+abs (a(nn-l,nn-2) ) 

x=0 .75*s 

y=x 

w= - 0 . 4 3 7 5 * s * * 2 
endif 
its=its+l 
do 15 m=nn -2,1, - 1 
z=a (m,m) 
r=x-z 
s=y-z 

p=(r*s-w) /a(m+l,m) +a (m,m-t-l) 

q-a (m+1 , m+1 ) - z - r - s 

r=a(m+2,m+l) 

s=abs (p) +abs (q) +abs ( r) 

p=p/s 

q=q/s 

r=r/s 

if (m.eq. 1) goto 4 

u=abs (a (m,m- 1) ) * (abs (q) +abs (r) ) 

v=abs (p) * (abs (a(m- l,m-l) ) +abs (z) +abs (a (m+l,m+l) ) ) 
if (u+v. eq. v) goto 4 
continue 
do 16 i=m+2,nn 
a(i, i-2) =0 . 

if (i.ne.m+2) a(i,i-3)=0. 
continue 
do 19 k=m,nn-l 
if (k.ne.m) then 
p=a (k, k- 1) 
q=a(k+l,k-l) 
r=0 . 
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if (k.ne.nn- 1) r=a (k+2 , k- 1) 
x=abs (p) +abs (q) +abs (r) 
if (x.ne.O . ) then 
p=p/x 
q=q/x 
r=r/x 
endif 
endif 

s=sign(sqrt (p**2+q**2+r**2) ,p) 
if (s .ne.O . ) then 
if (k.eq.m) then 

if (l.ne.m)a(k,k-l)=-a(k,k-l) 
else 

a(k,k-l)=-s*x 
endif 
p=p+s 
x=p/s 
y=q/s 
z=r/ s 

q=q/p 

r=r/p 

do 17 j=k,nn 

p=a (k, j ) +q*a (k+1, j ) 
if (k .ne.nn- 1) then 
p=p+r*a(k+2, j ) 
a (k+2, j ) =a (k+2, j ) -p*z 
endif 

a (k+1, j ) =a (k+1, j ) -p*y 
a (k, j ) =a (k, j ) -p*x 

17 continue 

do 18 i=l,min(nn,k+3) 
p=x*a (i, k) +y*a (i,k+l) 
if (k.ne.nn-1) then 
p=p+z*a(i,k+2) 
a ( i , k+2 ) =a ( i , k+2 ) - p* r 
endif 

a ( i, k+1) =a (i, k+1) -p*q 
a(i,k)=a(i,k) -p 

18 continue 
endif 

19 continue 
goto 2 

endif 
endif 
goto 1 
endif 
return 
EM) 

C (C) Copr. 1986-92 Numerical Recipes Software 41. 9 ;%{$!. 

SUBROUTINE sort2a (n, arr , brr) 
implicit real*8 (a-h,p-z) 
real*8 arr (n) , brr (n) 
real*8 tl (100) , t2 (100) 
do i - 1 , n 

tl(i) = arr(i) 
t2 (i) = brr (i) 
enddo 


12 



do i = 1 , n 
amin=arr ( i ) 
imin=i 

do j = i+l,n 

if ( (arr(j) .ge.O.and. (amin.lt.O.or.arr(j) .lt.amin) ) .or. 
. (amin. It. 0. and. arr(j) .lt.amin) ) then 

amin = arr(j) 
imin = j 
endif 
enddo 

if (imin.ne.i) then 
temp = arr(i) 
ar r ( i ) = ar r ( imin) 
arr(imin) = temp 
temp = brr(i) 
brr(i) = brr(imin) 
brr(imin) - temp 
endif 
enddo 

write(70,*) 'sort output:' 
do i = l,n 

write (70, ' <4el5.5) ' ) arr (i) ,brr (i) , tl (i) , t2 (i) 
enddo 
return 
END 

SUBROUTINE sort2 (n, arr , brr) 

INTEGER n, M, NSTACK 
REAL* 8 arr(n) ,brr(n) 

PARAMETER (M=7 , NSTACK=50) 

INTEGER i, ir, j , j stack, k, 1, istack (NSTACK) 

REAL* 8 a , b , temp 
j stack=0 
1=1 
ir=n 

1 if (ir-l.lt.M) then 

do 12 j=l+l,ir 
a=arr ( j ) 
b=brr ( j ) 

do 11 i=j -1,1, - 1 

if (arr(i) ,le.a)goto 2 
arr (i+1) =arr (i) 
brr(i+l)=brr(i) 

11 continue 
i=0 

2 arr (i+1) =a 
brr (i+1) =b 

12 continue 

if (jstack.eq.O) return 
ir=istack ( j stack) 
l=istack ( j stack-1) 
j s tack= j s tack - 2 
else 

k= (1+ir) /2 
temp=arr (k) 
arr (k) =arr (1+1) 
arr (1+1) =temp 
temp=brr (k) 
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brr (k) =brr ( 1+1 ) 
brr (l+l)=temp 

if (arr (1+1) .gt.arr (ir) ) then 
temp=arr ( 1+1 ) 
arr (1+1) =arr (ir) 
arr (ir)=temp 
temp=brr (1+1) 
brr (1+1) =brr (ir) 
brr (ir) =temp 
endif 

if (arr (1) .gt.arr (ir) ) then 
temp^arr (1) 
arr (l)=arr (ir) 
arr (ir) =temp 
temp=brr (1) 
brr (1) =brr (ir) 
brr (ir) =temp 
endif 

if (arr (1+1) .gt.arr (1) ) then 
temp=arr (1+1) 
arr (1+1) =arr (1) 
arr (1) -temp 
temp=brr (1+1) 
brr (1+1) =brr (1) 
brr(l)=temp 
endif 
i=l+l 
j=ir 
a=arr ( 1 ) 
b=brr (1) 

3 continue 

i=i+l 

if (arr (i) . It .a) goto 3 

4 continue 

j=j -1 

if (arr(j) .gt.a)goto 4 

if ( j . It . i) goto 5 

temp=arr (i) 

arr (i) =arr ( j ) 

arr ( j ) =temp 

temp=brr (i) 

brr(i)=brr(j) 

brr ( j ) =temp 

goto 3 

5 arr(l)=arr(j) 
arr ( j ) =a 
brr (1) =brr (j) 
brr ( j ) =b 

j stack= j stack+2 

if ( j stack. gt.NSTACK) pause ' NSTACK too small in sort2 ' 
if (ir- i+1 .ge. j - 1) then 
istack ( j stack) =ir 
istack (j stack- 1) =i 
ir=j - 1 
else 

istack (j stack) =j - 1 
istack (j stack- 1) =1 
l=i 
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endif 
endif 
goto 1 
END 

C (C) Copr. 1986-92 Numerical Recipes Software 41.9;%{$1. 

subroutine sysid (iter, d, v, n, dt , f f , stm, ldstm, iconv, q,nq, qcnaxq, 
. iflut,idone) 

implicit real*8 (a-h,p-z) 
parameter (nmax=50,nlag=l) 

real*8 cO (nmax, nmax) ,cl (nmax,nmax) ,dold(nmax) .void (nmax) 

real*8 a (nmax, nmax) , freq(nmax) , damp (nmax) 

real*8 freqold(nmax) ,dampold(nmax) .delta (nlag+1, nmax) 

real*8 d(n) ,v(n) ,p (nmax) ,aa,wr (nmax) ,wi (nmax) 

real*8 stm(ldstm, 2*n) , a2 (nmax, nmax) 

complex* 16 e 

LOGICAL sing 

integer ifirst 

save cO, cl, dold, void, ifirst 

data ifirst/1/ 

data tol/-l.d-6/ 

write (*,*) 'entered sysid. iter =',iter 

c initialize covariance matrices. 


c== : 

c 

c— — : 


c=== 


if (ifirst.eq.l) then 

write(*,*) 'in first loop.' 
ifirst = 0 
do i = 1 , nmax 
do j = l.nmax 
cO (i, j ) =0 
cl (i, j) =0 
enddo 
enddo 


update cO (the covariance matrix at dt=0) 


else 


write(*,*) 'updating cO ' 
do i = l,n 


do j = 1 , n 

cO (2*i - 1 , 2* j - 1) 
cO (2*i , 2* j -1) 

cO ( 2 * i - 1 , 2*j ) 

cO (2*i , 2 * j ) 

enddo 

cO (2*i-l, 2*n+l) = 
c0(2*i ,2*n+l) - 

cO (2*n+l, 2*i-l) = 
c0(2*n+l,2*i ) = 

enddo 


= ff*c0(2*i-l,2*j -1) +dold ( i ) *do!d(j) 
= f f *c0 (2*i , 2* j - 1) +vold (i) *dold ( j ) 
= f f *c0 (2*i - 1, 2* j ) +dold(i) *vold{ j ) 
= ff*c0(2*i ,2*j ) +vold(i) *vold( j ) 

cO (2*i - 1, 2*n+l) +dold(i) 
c0(2*i , 2*n+l) +vold(i) 
cO (2*n+l, 2*i - 1) +dold (i) 
c0(2*n+l,2*i )+vold(i) 


cO (2*n+l, 2*n+l) = cO (2 *n+l , 2 *n+l) +1 . 0 
do i = 2,2*n+l 
do j = 1 , i - 1 

cO (i, j ) = cO (j ,i) 
enddo 
enddo 
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c 


c= : 


c : 

c 

c- 


C-- 

c 

c= : 


update cl (the transpose of the covariance matrix at dt=l) 


write(*,*) 'updating cl' 
do i = l,n 


do j = l,n 

cl (2*i-l,2*j -1) 
cl(2*i / 2*j -1) 
Cl (2*i-l,2*j ) 

cl (2*i ,2*i ) 

enddo 


f f *cl (2*i-l,2*j-l) +dold (i) *d( j ) 
ff*cl(2*i , 2*j -1) +vold(i) *d ( j ) 
ff *cl (2*i-l, 2*j ) +dold(i) *v(j ) 

ff*cl(2*i , 2*3 )+vold(i) *v(j) 


cl (2*i -1, 2*n+l) 
cl (2*i ,2*n+l) 

cl (2*n+l, 2*i-l) 
cl(2*n+l,2*i ) 

enddo 

cl (2*n+l, 2*n+l) = 


= cl (2*i -1, 2*n+l) +dold (i) 
= cl{2*i , 2*n+l) +vold(i) 
= cl (2*n+l, 2*i-l) +d (i) 

= cl(2*n+l,2*i )+v(i) 

cl(2*n+l / 2*n+l)+1.0 


endif 


store the displacement and velocity information. 


write{19, *) 'd,v: ' 
do i = l,n 

write{19, ' (2el2.5) '} dold(i),d(i) 
write (19, ' (2el2 .5) ’ ) vold(i),v(i) 
enddo 

do i = l,n 

dold(i) = d(i) 
vold(i) = v(i) 
enddo 

write (19 , *) 'cO : 1 
do i = l,2*n+l 

write (19 , ' (255el2.5) ') (cO (i, j ) , j=l, 2*n+l) 
enddo 

write(19, *) 'cl: ' 
do i = l,2*n+l 

write (19 , ' (255el2.5) ') (cl(i,j) , j=l,2*n+l) 
enddo 


perform the backsubstitution to get A (if cO is non-singular) . 


call choldc(c0,2*n+l,ninax,p,sing) 
write (*, *) 'qO : 1 
do i = l,2*n 

write (*, ' (255el5 .8) ') (0, j=l,i-l) ,p(i) , (cO (j , i) , j=i+l, 2*n) 

enddo 

if (sing) goto 900 

write(*,*) 'backsubstituting. ' 

call cholsl (cO, 2*n+l, nmax, 2*n, p, cl, a) 

wri te ( * , * ) ' a : ' 

do i = l,2*n+l 

write (*, ’ (255el5.8) ') (a(i,j) , j=l, 2*n+l) 

enddo 

do i = l,2*n 

do j = i+l,2*n+l 
aa = a (i, j ) 
a (i, j ) = a (j , i) 
a(j ,i) = a a 
enddo 
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enddo 

write {*, *) 'a: 1 
do i = l,2*n+l 

write (*, ' (255el5 .8) ’ ) (a (i, j ) , j=l, 2*n+l) 
enddo 

c find the frequencies and damping ratios. 

write(*,*) 'computing frequencies.' 
call elmhes(a,2*n,nmax) 
call hqr (a, 2*n,nmax,wr,wi) 
call sort2a(2*n,wi,wr) 
do i = 2, nlag 
do j = l,n 

delta (i, j ) =delta(i-l, j) 
enddo 
enddo 

do i - l,n 

e = zlog(dcmplx(wr (i) ,wi (i) ) ) /dt 
freq(i) = abs (imag (e) ) /2/3 . 14159 
damp(i) = real (e) /abs (e) 

delta (1, i) = sqrt ( { (f req (i) -f reqold (i) ) /f req (i) ) **2+ 
( (damp(i) -dampold(i) )/damp(i) ) **2) 

enddo 

do j = 1 , n 

delta (nlag+1 , j ) = 0 
do i = l,nlag 

delta (nlag+1 , j ) = delta (nlag+1 , j ) +delta (i , j ) 
enddo 
enddo 

write (16, ' (255el5 . 5) ' ) (wr (i) ,wi (i) , i=n+l,2*n) 
c check convergence of the frequencies & damping ratios, 
dd = 0 

do i = n+l,2*n 

dd = dd+delta (nlag+1, i) 
enddo 

dd = dd/nlag/n 

write (iconv, ' (255el3 .5) ' ) (freq(i) ,damp(i) , i=l,n) , 
(delta(nlag+l, i) , i=l, n),dd 
do i = n+l,2*n 

freqold(i) = freq(i) 
dampold(i) = damp(i) 
enddo 

if ( (dd.ge. 0 . OdO .and. dd. It. tol) .or . idone.ne.O) then 

c reconstruct the A matrix. 

call cholsl (cO, 2*n+l,nmax, 2*n+l, p, cl, a) 
do i = l,2*n-l 
do j = i+1, 2*n 
aa = a (i, j ) 
a ( i , j ) = a(j ,i) 
a ( j , i) = aa 
enddo 
enddo 

write (*, *) 'a: ' 
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do i = l,2*n 

write (*, ' (255el5 .8) ') (a (i, j ) , j=l, 2*n) 
enddo 

c do a flutter analysis to get q-g and q-f plots. 

do iq= G,nq 

qq = iq* (q*qmaxq/nq) 
do j = l,2*n 
do k = l,2*n 

a2(j,k) = stm(j ,k) + (a(j ,k) -stm(j ,k) ) *qq/q 
enddo 
enddo 

write (*,*) 'a2: ' 
do i = l,2*n 

write (*, ' (255el5 .5) ') (a2 (i, j ) , j=l, 2*n) 
enddo 

call elmhes (a2, 2*n,nmax) 
call hqr (a2, 2*n,nmax,wr,wi) 
do i = l,2*n 

write(*,*) wr(i),wi(i) 
e = zlog(dcmplx(wr (i) ,wi (i) ) ) /dt 
wi (i) = imag(e)/2/3. 14159 
wr(i) = real (e) /abs (e) *2 
enddo 

write(70,*) ' iq = 1 , iq 
call sort2a (2*n, wi, wr) 
do i = l,n 

f req (i) = wi (i) 
damp ( i ) = wr ( i ) 
enddo 

write(20, ' (255el5.5) ') qq*144, (freq(i) ,damp(i) ,i=l,n) 
if (iq.gt.l) then 
do i = l,n 

if (damp (i) *dampold (i) . It .0) then 

qqq = qq- (q*qmaxq/nq) *damp (i) / (damp (i) -dampold (i) ) 
f = freq(i) - (freq(i) -freqold(i) ) *damp(i) / 

. (dait^) (i) -dampold (i) ) 

write (if lut, ' (26hflutter crossing for mode ,i3,6h at q-, 
f8.4,7h and f=,f8.4)') i-n,qqq,f 
endif 

dampold (i) = damp(i) 
freqold(i) = freq(i) 
enddo 
endif 
enddo 
idone = 1 
endif 

900 continue 
return 
end 
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LIMITED EXCLUSIVE RIGHTS NOTICE 
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Figure 1 





Time History of Generalized Displacements for 
Flexible Semispan Model (FSM) 

HSCT Aerodynamics, Long Beach 



CFL3D.AE-BA, Euler, = 0.80, a = 2.0°, q = 250 psf, 65x241x49 C~0 Grid 
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Q-F Plot Computed From Linear Flutter Analysis 
Mach 0.80, a=2.0 deg, q0=250 PSF. 



Spec. Doc. Collection 


Q-G Plot Computed From Nonlinear Flutter Simulation 
Mach 0.80, a=2.0 deg, q0=250 PSF. 


Exact Solutionsj at 
Q=0 and Q=250 







